#include "StdAfx.h"
#include "WaveletTransform.h"
#include <cmath>
#define Haar_Wavelet 1
#define Daub_Wavelet 2

CWaveletTransform::CWaveletTransform(double** input, int inputcolumn, int inputrow)
{
	m_InputData = input;
	m_InputColumn = inputcolumn;
	m_InputRow = inputrow;
}

CWaveletTransform::~CWaveletTransform(void)
{
}

//Set the parameters of the wavelet transform and create memory for each level
bool CWaveletTransform::SetParameters(int DecomposeLevel, int TypeOfWavelet)
{
	//Validate the level of decomposition
	if(DecomposeLevel<=0)
	{
		::AfxMessageBox(CString("Man, don't joke!!! Make the decomposing level >= 1"), MB_OK, 0);
		return 0;
	}
	if((m_InputColumn/(int)pow(2.0, DecomposeLevel))<1||(m_InputRow/(int)pow(2.0, DecomposeLevel))<2)
	{
		::AfxMessageBox(CString("Decompose level too much, select a smaller number"), MB_OK, 0);
		return 0;
	}

	//Set the parameters of the wavelet transform
	m_DecompositionLevel = DecomposeLevel;
	m_TypeOfWavelet = TypeOfWavelet;

	//Calculate the column and row at each level
	m_NewColumn = new int[m_DecompositionLevel];
	m_NewRow = new int[m_DecompositionLevel];
	m_NewColumn[0] = m_InputColumn/2;
	m_NewRow[0] = m_InputRow/2;
	for(int i=1; i<m_DecompositionLevel; i++)
	{
		m_NewColumn[i] = m_NewColumn[i-1]/2;
		m_NewRow[i] = m_NewRow[i-1]/2;
	}

	//Create memory for each level
	Approximation = new double**[m_DecompositionLevel];
	Horizontal = new double**[m_DecompositionLevel];
	Vertical = new double**[m_DecompositionLevel];
	Diagonal = new double**[m_DecompositionLevel];
	for(int k=0; k<m_DecompositionLevel; k++)
	{
		Approximation[k] = new double*[m_NewRow[k]];
		Horizontal[k] = new double*[m_NewRow[k]];
		Vertical[k] = new double*[m_NewRow[k]];
		Diagonal[k] = new double*[m_NewRow[k]];
		for(int i=0; i<m_NewRow[k]; i++)
		{
			Approximation[k][i] = new double[m_NewColumn[k]];
			Horizontal[k][i] = new double[m_NewColumn[k]];
			Vertical[k][i] = new double[m_NewColumn[k]];
			Diagonal[k][i] = new double[m_NewColumn[k]];
		}
	}
	return 1;
}

//Get the approximation at certain level
double** CWaveletTransform::GetApproximation(int ApproximationLevel, int *DataColumn, int *DataRow)
{
	//Validate the levels
	if(ApproximationLevel>m_DecompositionLevel||ApproximationLevel<=0)
	{
		::AfxMessageBox(CString("Illegal approximation level, first entry is returned"), MB_OK, 0);
		*DataColumn = 0;
		*DataRow = 0;
		return Approximation[0];
	}
	*DataColumn = m_NewColumn[ApproximationLevel-1];
	*DataRow = m_NewRow[ApproximationLevel-1];
	return Approximation[ApproximationLevel-1];
}

//Get the horizontal decomposition
double** CWaveletTransform::GetHorizontal(int HorizontalLevel, int *DataColumn, int *DataRow)
{
	//Validate the levels
	if(HorizontalLevel>m_DecompositionLevel||HorizontalLevel<=0)
	{
		::AfxMessageBox(CString("Illegal horizontal level, first entry is returned"), MB_OK, 0);
		*DataColumn = 0;
		*DataRow = 0;
		return Horizontal[0];
	}
	*DataColumn = m_NewColumn[HorizontalLevel-1];
	*DataRow = m_NewRow[HorizontalLevel-1];
	return Horizontal[HorizontalLevel-1];
}

//Get the vertical decomposition
double** CWaveletTransform::GetVertical(int VerticalLevel, int *DataColumn, int *DataRow)
{
	//Validate the levels
	if(VerticalLevel>m_DecompositionLevel||VerticalLevel<=0)
	{
		::AfxMessageBox(CString("Illegal vertical level, first entry is returned"), MB_OK, 0);
		*DataColumn = 0;
		*DataRow = 0;
		return Vertical[0];
	}
	*DataColumn = m_NewColumn[VerticalLevel-1];
	*DataRow = m_NewRow[VerticalLevel-1];
	return Vertical[VerticalLevel-1];
}

//Get the diagonal decomposition
double** CWaveletTransform::GetDiagonal(int DiagonalLevel, int *DataColumn, int *DataRow)
{
	//Validate the levels
	if(DiagonalLevel>m_DecompositionLevel||DiagonalLevel<=0)
	{
		::AfxMessageBox(CString("Illegal diagonal level, first entry is returned"), MB_OK, 0);
		*DataColumn = 0;
		*DataRow = 0;
		return Diagonal[0];
	}
	*DataColumn = m_NewColumn[DiagonalLevel-1];
	*DataRow = m_NewRow[DiagonalLevel-1];
	return Diagonal[DiagonalLevel-1];
}

//Perform the wavelet transform
void CWaveletTransform::PerformTransform()
{
	//Wavelet transform of the first level
	int decomposelevel = 0;
	CWaveletTransform::WaveletTransform(m_InputData, m_InputColumn, m_InputRow, Approximation[0],
		Horizontal[0], Vertical[0], Diagonal[0], &m_NewColumn[0], &m_NewRow[0]);
	decomposelevel ++;
	//Wavelet transform at different levels
	for(;decomposelevel<m_DecompositionLevel; decomposelevel++)
	{
		CWaveletTransform::WaveletTransform(Approximation[decomposelevel-1],
			m_NewColumn[decomposelevel-1], m_NewRow[decomposelevel-1],
			Approximation[decomposelevel], Horizontal[decomposelevel],
			Vertical[decomposelevel], Diagonal[decomposelevel],
			&m_NewColumn[decomposelevel], &m_NewRow[decomposelevel]);
	}
	return;
}

//Perform a wavelet transform based on the input
bool CWaveletTransform::WaveletTransform(double **input, int inputcolumn, int inputrow,
										 double **approximation, double **horizontal,
										 double **vertical, double **diagonal,
										 int *newdatacolumn,int *newdatarow)
{
	switch(m_TypeOfWavelet)
	{
	case Haar_Wavelet:
		break;
	case Daub_Wavelet:
		{
			//Define the filter
			double LowFilter[4] = {(1.0+sqrt(3.0))/(4.0*sqrt(2.0)), (3.0+sqrt(3.0))/(4.0*sqrt(2.0)), 
				(3.0-sqrt(3.0))/(4.0*sqrt(2.0)), (1.0-sqrt(3.0))/(4.0*sqrt(2.0))};
			double HighFilter[4] = {LowFilter[3], -LowFilter[2], LowFilter[1], -LowFilter[0]};
			//Create memory for the low and high pass filter
			double** LowPassFilterMatrix = new double*[inputrow];
			double** HighPassFilterMatrix = new double*[inputrow];
			for(int i=0; i<inputrow; i++)
			{
				LowPassFilterMatrix[i] = new double[inputcolumn/2];
				HighPassFilterMatrix[i] = new double[inputcolumn/2];
			}
			//Low pass filter by rows and downsampling by columns
			CWaveletTransform::WaveletFilterRow(input, inputcolumn, inputrow, 
				LowPassFilterMatrix, &LowFilter[0], 4);
			//High pass filter by rows and downsampling by columns
			CWaveletTransform::WaveletFilterRow(input, inputcolumn, inputrow, 
				HighPassFilterMatrix, &HighFilter[0], 4);
			//Low pass filter by columns and downsampling by rows of the low pass filtered matrix
			CWaveletTransform::WaveletFilterColumn(LowPassFilterMatrix,	inputcolumn/2, inputrow,
				approximation, &LowFilter[0], 4);
			//High pass filter by columns and downsampling by rows of the low pass filter matrix
			CWaveletTransform::WaveletFilterColumn(LowPassFilterMatrix,	inputcolumn/2, inputrow,
				horizontal, &HighFilter[0], 4);
			//Low pass filter by columns and downsampleing by rows of the high pass filter matrix
			CWaveletTransform::WaveletFilterColumn(HighPassFilterMatrix, inputcolumn/2, inputrow,
				vertical, &LowFilter[0], 4);
			//High pass filter by columns and downsampling by rows of the high pass filter matrix
			CWaveletTransform::WaveletFilterColumn(HighPassFilterMatrix, inputcolumn/2, inputrow,
				diagonal, &HighFilter[0], 4);
			//Delete the old matrix
			bf.DeleteMatrix(LowPassFilterMatrix, inputrow);
			bf.DeleteMatrix(HighPassFilterMatrix, inputrow);
		}
		break;
	default:
		break;
	}
	return 1;
}

//Filter a matrix with the wavelet filter in columns and downsampling by rows
void CWaveletTransform::WaveletFilterColumn(double **input, int inputcolumn, int inputrow,
												double **destination, double *filter, int filtersize)
{
	//Define a contrainer for filter data
	double* filterdata = new double[filtersize];
	//Filter by columns
	for(int i=0; i<inputrow/2; i++)
	{
		for(int j=0; j<inputcolumn; j++)
		{
			destination[i][j] = 0;
			//Filtering
			for(int k=0; k<filtersize; k++)
			{
				//Find a proper value for the current filter data
				if((i*2-filtersize/2+k)<0)
					filterdata[k] = input[inputrow+(i*2-filtersize/2+k)][j];
				else
					if((i*2-filtersize/2+k)>=inputrow)
						filterdata[k] = input[(i*2-filtersize/2+k)-inputrow][j];
					else
						filterdata[k] = input[i*2-filtersize/2+k][j];
				destination[i][j] += filter[k]*filterdata[k];
			}
		}
	}
	//Clean memory and return;
	free(filterdata);
	return;
}

//Wavelet filter by rows and downsampling by columns
void CWaveletTransform::WaveletFilterRow(double **input, int inputcolumn, int inputrow, 
											 double **destination, double *filter, int filtersize)
{
	//Define a container for the filter data
	double* filterdata = new double[filtersize];
	//Filtering by rows
	for(int i=0; i<inputrow; i++)
	{
		for(int j=0; j<inputcolumn/2; j++)
		{
			destination[i][j] = 0;
			//Filtering
			for(int k=0; k<filtersize; k++)
			{
				//Determine the value for the current filterdata
				if((j*2-filtersize/2+k)<0)
					filterdata[k] = input[i][inputcolumn+(j*2-filtersize/2+k)];
				else
					if((j*2-filtersize/2+k)>=inputcolumn)
						filterdata[k] = input[i][(j*2-filtersize/2+k)-inputcolumn];
					else
						filterdata[k] = input[i][j*2-filtersize/2+k];
				destination[i][j] += filter[k]*filterdata[k];
			}
		}
	}
	//Clean memory and return;
	free(filterdata);
	return;
}

//Destroy the wavelet transform
void CWaveletTransform::WaveletDestroy(void)
{
	//Clean the data in each level
	for(int i=0; i<m_DecompositionLevel; i++)
	{
		for(int j=0; j<m_NewRow[i]; j++)
		{
			free(Approximation[i][j]);
			free(Horizontal[i][j]);
			free(Vertical[i][j]);
			free(Diagonal[i][j]);
		}
		free(Approximation[i]);
		free(Horizontal[i]);
		free(Vertical[i]);
		free(Diagonal[i]);
	}
	//Clean the number of columns and rows at each level
	free(m_NewColumn);
	free(m_NewRow);
}

//Get vertical decompositions
double*** CWaveletTransform::GetVertical(void)
{
	return Vertical;
}

//Get the horizontal decompositions
double*** CWaveletTransform::GetHorizontal(void)
{
	return Horizontal;
}

//Get the diagonal decompositions
double*** CWaveletTransform::GetDiagonal(void)
{
	return Diagonal;
}

//Get the approximations
double*** CWaveletTransform::GetApproximation(void)
{
	return Approximation;
}

//Get the columns
int* CWaveletTransform::GetColumns(void)
{
	return m_NewColumn;
}

//Get the rows
int* CWaveletTransform::GetRows(void)
{
	return m_NewRow;
}
